{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimizing DiffEq Code\n",
    "\n",
    "In this notebook we will walk through some of the main tools for optimizing your code in order to efficiently solve DifferentialEquations.jl. User-side optimizations are important because, for sufficiently difficult problems, most of the time will be spent inside of your `f` function, the function you are trying to solve. \"Efficient\" integrators are those that reduce the required number of `f` calls to hit the error tolerance. The main ideas for optimizing your DiffEq code, or any Julia function, are the following:\n",
    "\n",
    "- Make it non-allocating\n",
    "- Use StaticArrays for small arrays\n",
    "- Use broadcast fusion\n",
    "- Make it type-stable\n",
    "- Reduce redundant calculations\n",
    "- Make use of BLAS calls\n",
    "- Optimize algorithm choice\n",
    "\n",
    "We'll discuss these strategies in the context of small and large systems. Let's start with small systems.\n",
    "\n",
    "## Optimizing Small Systems (<100 DEs)\n",
    "\n",
    "Let's take the classic Lorenz system from before. Let's start by naively writing the system in its out-of-place form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lorenz (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function lorenz(u,p,t)\n",
    " dx = 10.0*(u[2]-u[1])\n",
    " dy = u[1]*(28.0-u[3]) - u[2]\n",
    " dz = u[1]*u[2] - (8/3)*u[3]\n",
    " [dx,dy,dz]\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, `lorenz` returns an object, `[dx,dy,dz]`, which is created within the body of `lorenz`.\n",
    "\n",
    "This is a common code pattern from high-level languages like MATLAB, SciPy, or R's deSolve. However, the issue with this form is that it allocates a vector, `[dx,dy,dz]`, at each step. Let's benchmark the solution process with this choice of function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  11.11 MiB\n",
       "  allocs estimate:  124161\n",
       "  --------------\n",
       "  minimum time:     4.649 ms (0.00% GC)\n",
       "  median time:      8.045 ms (41.31% GC)\n",
       "  mean time:        7.412 ms (28.65% GC)\n",
       "  maximum time:     26.596 ms (67.17% GC)\n",
       "  --------------\n",
       "  samples:          674\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using DifferentialEquations, BenchmarkTools\n",
    "u0 = [1.0;0.0;0.0]\n",
    "tspan = (0.0,100.0)\n",
    "prob = ODEProblem(lorenz,u0,tspan)\n",
    "@benchmark solve(prob,Tsit5())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The BenchmarkTools package's `@benchmark` runs the code multiple times to get an accurate measurement. The minimum time is the time it takes when your OS and other background processes aren't getting in the way. Notice that in this case it takes about 5ms to solve and allocates around 11.11 MiB. However, if we were to use this inside of a real user code we'd see a lot of time spent doing garbage collection (GC) to clean up all of the arrays we made. Even if we turn off saving we have these allocations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  9.77 MiB\n",
       "  allocs estimate:  111505\n",
       "  --------------\n",
       "  minimum time:     4.014 ms (0.00% GC)\n",
       "  median time:      5.043 ms (0.00% GC)\n",
       "  mean time:        5.667 ms (18.53% GC)\n",
       "  maximum time:     107.954 ms (1.88% GC)\n",
       "  --------------\n",
       "  samples:          881\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark solve(prob,Tsit5(),save_everystep=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem of course is that arrays are created every time our derivative function is called. This function is called multiple times per step and is thus the main source of memory usage. To fix this, we can use the in-place form to ***make our code non-allocating***:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lorenz! (generic function with 1 method)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function lorenz!(du,u,p,t)\n",
    " du[1] = 10.0*(u[2]-u[1])\n",
    " du[2] = u[1]*(28.0-u[3]) - u[2]\n",
    " du[3] = u[1]*u[2] - (8/3)*u[3]\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, instead of creating an array each time, we utilized the cache array `du`. When the inplace form is used, DifferentialEquations.jl takes a different internal route that minimizes the internal allocations as well. When we benchmark this function, we will see quite a difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  1.33 MiB\n",
       "  allocs estimate:  12612\n",
       "  --------------\n",
       "  minimum time:     1.230 ms (0.00% GC)\n",
       "  median time:      1.406 ms (0.00% GC)\n",
       "  mean time:        1.670 ms (10.92% GC)\n",
       "  maximum time:     9.507 ms (78.55% GC)\n",
       "  --------------\n",
       "  samples:          2980\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u0 = [1.0;0.0;0.0]\n",
    "tspan = (0.0,100.0)\n",
    "prob = ODEProblem(lorenz!,u0,tspan)\n",
    "@benchmark solve(prob,Tsit5())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  8.34 KiB\n",
       "  allocs estimate:  126\n",
       "  --------------\n",
       "  minimum time:     711.216 μs (0.00% GC)\n",
       "  median time:      820.568 μs (0.00% GC)\n",
       "  mean time:        866.374 μs (0.17% GC)\n",
       "  maximum time:     5.141 ms (82.17% GC)\n",
       "  --------------\n",
       "  samples:          5724\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark solve(prob,Tsit5(),save_everystep=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a 4x time difference just from that change! Notice there are still some allocations and this is due to the construction of the integration cache. But this doesn't scale with the problem size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  8.34 KiB\n",
       "  allocs estimate:  126\n",
       "  --------------\n",
       "  minimum time:     3.964 ms (0.00% GC)\n",
       "  median time:      4.252 ms (0.00% GC)\n",
       "  mean time:        4.504 ms (0.00% GC)\n",
       "  maximum time:     11.964 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          1108\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tspan = (0.0,500.0) # 5x longer than before\n",
    "prob = ODEProblem(lorenz!,u0,tspan)\n",
    "@benchmark solve(prob,Tsit5(),save_everystep=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "since that's all just setup allocations. \n",
    "\n",
    "#### But if the system is small we can optimize even more.\n",
    "\n",
    "Allocations are only expensive if they are \"heap allocations\". For a more in-depth definition of heap allocations, [there are a lot of sources online](http://net-informations.com/faq/net/stack-heap.htm). But a good working definition is that heap allocations are variable-sized slabs of memory which have to be pointed to, and this pointer indirection costs time. Additionally, the heap has to be managed and the garbage controllers has to actively keep track of what's on the heap.\n",
    "\n",
    "However, there's an alternative to heap allocations, known as stack allocations. The stack is statically-sized (known at compile time) and thus its accesses are quick. Additionally, the exact block of memory is known in advance by the compiler, and thus re-using the memory is cheap. This means that allocating on the stack has essentially no cost!\n",
    "\n",
    "Arrays have to be heap allocated because their size (and thus the amount of memory they take up) is determined at runtime. But there are structures in Julia which are stack-allocated. `struct`s for example are stack-allocated \"value-type\"s. `Tuple`s are a stack-allocated collection. The most useful data structure for DiffEq though is the `StaticArray` from the package [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl). These arrays have their length determined at compile-time. They are created using macros attached to normal array expressions, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element SVector{3,Float64}:\n",
       " 2.0\n",
       " 3.0\n",
       " 5.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using StaticArrays\n",
    "A = @SVector [2.0,3.0,5.0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the `3` after `SVector` gives the size of the `SVector`. It cannot be changed. Additionally, `SVector`s are immutable, so we have to create a new `SVector` to change values. But remember, we don't have to worry about allocations because this data structure is stack-allocated. `SArray`s have a lot of extra optimizations as well: they have fast matrix multiplication, fast QR factorizations, etc. which directly make use of the information about the size of the array. Thus, when possible they should be used.\n",
    "\n",
    "Unfortunately static arrays can only be used for sufficiently small arrays. After a certain size, they are forced to heap allocate after some instructions and their compile time balloons. Thus static arrays shouldn't be used if your system has more than 100 variables. Additionally, only the native Julia algorithms can fully utilize static arrays.\n",
    "\n",
    "Let's ***optimize `lorenz` using static arrays***. Note that in this case, we want to use the out-of-place allocating form, but this time we want to output a static array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lorenz_static (generic function with 1 method)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function lorenz_static(u,p,t)\n",
    " dx = 10.0*(u[2]-u[1])\n",
    " dy = u[1]*(28.0-u[3]) - u[2]\n",
    " dz = u[1]*u[2] - (8/3)*u[3]\n",
    " @SVector [dx,dy,dz]\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make the solver internally use static arrays, we simply give it a static array as the initial condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  408.94 KiB\n",
       "  allocs estimate:  2249\n",
       "  --------------\n",
       "  minimum time:     740.817 μs (0.00% GC)\n",
       "  median time:      823.518 μs (0.00% GC)\n",
       "  mean time:        887.751 μs (3.49% GC)\n",
       "  maximum time:     5.072 ms (82.55% GC)\n",
       "  --------------\n",
       "  samples:          5582\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u0 = @SVector [1.0,0.0,0.0]\n",
    "tspan = (0.0,100.0)\n",
    "prob = ODEProblem(lorenz_static,u0,tspan)\n",
    "@benchmark solve(prob,Tsit5())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  8.95 KiB\n",
       "  allocs estimate:  121\n",
       "  --------------\n",
       "  minimum time:     576.313 μs (0.00% GC)\n",
       "  median time:      687.015 μs (0.00% GC)\n",
       "  mean time:        755.745 μs (0.17% GC)\n",
       "  maximum time:     100.796 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          6551\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark solve(prob,Tsit5(),save_everystep=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that's pretty much all there is to it. With static arrays you don't have to worry about allocating, so use operations like `*` and don't worry about fusing operations (discussed in the next section). Do \"the vectorized code\" of R/MATLAB/Python and your code in this case will be fast, or directly use the numbers/values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 1\n",
    "\n",
    "Implement the out-of-place array, in-place array, and out-of-place static array forms for the [Henon-Heiles System](https://en.wikipedia.org/wiki/H%C3%A9non%E2%80%93Heiles_system) and time the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimizing Large Systems\n",
    "\n",
    "### Interlude: Managing Allocations with Broadcast Fusion\n",
    "\n",
    "When your system is sufficiently large, or you have to make use of a non-native Julia algorithm, you have to make use of `Array`s. In order to use arrays in the most efficient manner, you need to be careful about temporary allocations. Vectorized calculations naturally have plenty of temporary array allocations. This is because a vectorized calculation outputs a vector. Thus:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  15.26 MiB\n",
       "  allocs estimate:  4\n",
       "  --------------\n",
       "  minimum time:     5.923 ms (30.94% GC)\n",
       "  median time:      8.066 ms (19.37% GC)\n",
       "  mean time:        8.961 ms (13.41% GC)\n",
       "  maximum time:     109.321 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          557\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(1000,1000); B = rand(1000,1000); C = rand(1000,1000)\n",
    "test(A,B,C) = A + B + C\n",
    "@benchmark test(A,B,C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That expression `A + B + C` creates 2 arrays. It first creates one for the output of `A + B`, then uses that result array to `+ C` to get the final result. 2 arrays! We don't want that! The first thing to do to fix this is to use broadcast fusion. [Broadcast fusion](https://julialang.org/blog/2017/01/moredots) puts expressions together. For example, instead of doing the `+` operations separately, if we were to add them all at the same time, then we would only have a single array that's created. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  7.63 MiB\n",
       "  allocs estimate:  6\n",
       "  --------------\n",
       "  minimum time:     3.281 ms (0.00% GC)\n",
       "  median time:      4.436 ms (0.00% GC)\n",
       "  mean time:        4.810 ms (10.98% GC)\n",
       "  maximum time:     8.350 ms (18.21% GC)\n",
       "  --------------\n",
       "  samples:          1037\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test2(A,B,C) = map((a,b,c)->a+b+c,A,B,C)\n",
    "@benchmark test2(A,B,C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Puts the whole expression into a single function call, and thus only one array is required to store output. This is the same as writing the loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  7.63 MiB\n",
       "  allocs estimate:  2\n",
       "  --------------\n",
       "  minimum time:     2.165 ms (0.00% GC)\n",
       "  median time:      2.493 ms (0.00% GC)\n",
       "  mean time:        3.196 ms (16.54% GC)\n",
       "  maximum time:     104.956 ms (1.50% GC)\n",
       "  --------------\n",
       "  samples:          1555\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function test3(A,B,C)\n",
    "    D = similar(A)\n",
    "    @inbounds for i in eachindex(A)\n",
    "        D[i] = A[i] + B[i] + C[i]\n",
    "    end\n",
    "    D\n",
    "end\n",
    "@benchmark test3(A,B,C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, Julia's broadcast is syntactic sugar for this. If multiple expressions have a `.`, then it will put those vectorized operations together. Thus:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  7.63 MiB\n",
       "  allocs estimate:  57\n",
       "  --------------\n",
       "  minimum time:     2.419 ms (0.00% GC)\n",
       "  median time:      2.761 ms (0.00% GC)\n",
       "  mean time:        3.401 ms (15.93% GC)\n",
       "  maximum time:     102.888 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          1463\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test4(A,B,C) = A .+ B .+ C\n",
    "@benchmark test4(A,B,C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is a version with only 1 array created (the output). Note that `.`s can be used with function calls as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100-element Array{Float64,1}:\n",
       " 1.14031 \n",
       " 0.437576\n",
       " 1.18087 \n",
       " 1.07615 \n",
       " 1.20897 \n",
       " 0.821996\n",
       " 1.28333 \n",
       " 0.233329\n",
       " 1.24629 \n",
       " 0.928158\n",
       " 1.29115 \n",
       " 0.699507\n",
       " 1.51438 \n",
       " ⋮       \n",
       " 0.42712 \n",
       " 1.50921 \n",
       " 0.515603\n",
       " 1.05544 \n",
       " 0.307876\n",
       " 1.36608 \n",
       " 0.701174\n",
       " 0.493429\n",
       " 0.302578\n",
       " 0.870267\n",
       " 1.5601  \n",
       " 0.99526 "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sin.(A) .+ sin.(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, the `@.` macro applys a dot to every operator: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  7.63 MiB\n",
       "  allocs estimate:  57\n",
       "  --------------\n",
       "  minimum time:     2.513 ms (0.00% GC)\n",
       "  median time:      2.894 ms (0.00% GC)\n",
       "  mean time:        3.475 ms (16.39% GC)\n",
       "  maximum time:     15.745 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          1432\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test5(A,B,C) = @. A + B + C #only one array allocated\n",
    "@benchmark test5(A,B,C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using these tools we can get rid of our intermediate array allocations for many vectorized function calls. But we are still allocating the output array. To get rid of that allocation, we can instead use mutation. Mutating broadcast is done via `.=`. For example, if we pre-allocate the output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = zeros(1000,1000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can keep re-using this cache for subsequent calculations. The mutating broadcasting form is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  160 bytes\n",
       "  allocs estimate:  4\n",
       "  --------------\n",
       "  minimum time:     1.964 ms (0.00% GC)\n",
       "  median time:      2.208 ms (0.00% GC)\n",
       "  mean time:        2.267 ms (0.00% GC)\n",
       "  maximum time:     6.596 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          2188\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test6!(D,A,B,C) = D .= A .+ B .+ C #only one array allocated\n",
    "@benchmark test6!(D,A,B,C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we use `@.` before the `=`, then it will turn it into `.=`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  160 bytes\n",
       "  allocs estimate:  4\n",
       "  --------------\n",
       "  minimum time:     1.972 ms (0.00% GC)\n",
       "  median time:      2.197 ms (0.00% GC)\n",
       "  mean time:        2.266 ms (0.00% GC)\n",
       "  maximum time:     4.096 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          2189\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test7!(D,A,B,C) = @. D = A + B + C #only one array allocated\n",
    "@benchmark test7!(D,A,B,C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that in this case, there is no \"output\", and instead the values inside of `D` are what are changed (like with the DiffEq inplace function). Many Julia functions have a mutating form which is denoted with a `!`. For example, the mutating form of the `map` is `map!`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  32 bytes\n",
       "  allocs estimate:  1\n",
       "  --------------\n",
       "  minimum time:     2.011 ms (0.00% GC)\n",
       "  median time:      2.374 ms (0.00% GC)\n",
       "  mean time:        2.615 ms (0.00% GC)\n",
       "  maximum time:     102.888 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          1902\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test8!(D,A,B,C) = map!((a,b,c)->a+b+c,D,A,B,C)\n",
    "@benchmark test8!(D,A,B,C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some operations require using an alternate mutating form in order to be fast. For example, matrix multiplication via `*` allocates a temporary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  7.63 MiB\n",
       "  allocs estimate:  2\n",
       "  --------------\n",
       "  minimum time:     201.438 ms (0.00% GC)\n",
       "  median time:      297.936 ms (0.00% GC)\n",
       "  mean time:        289.920 ms (0.18% GC)\n",
       "  maximum time:     547.904 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          18\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark A*B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead, we can use the mutating form `A_mul_B!` into a cache array to avoid allocating the output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  0 bytes\n",
       "  allocs estimate:  0\n",
       "  --------------\n",
       "  minimum time:     200.741 ms (0.00% GC)\n",
       "  median time:      202.884 ms (0.00% GC)\n",
       "  mean time:        253.591 ms (0.00% GC)\n",
       "  maximum time:     398.854 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          21\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark A_mul_B!(D,A,B) # same as D = A * B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For repeated calculations this reduced allocation can stop GC cycles and thus lead to more efficient code. Additionally, ***we can fuse together higher level linear algebra operations using BLAS***. The package [SugarBLAS.jl](https://github.com/lopezm94/SugarBLAS.jl) makes it easy to write higher level operations like `alpha*B*A + beta*C` as mutating BLAS calls. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example Optimization: Gierer-Meinhardt Reaction-Diffusion PDE Discretization\n",
    "\n",
    "Let's optimize the solution of a Reaction-Diffusion PDE's discretization. In its discretized form, this is the ODE:\n",
    "\n",
    "$$ du = D_1 (A_y u + u A_x) + \\frac{au^2}{v} + \\bar{u} - \\alpha u $$\n",
    "$$ dv = D_2 (A_y v + v A_x) + a u^2 + \\beta v $$\n",
    "\n",
    "where $u$, $v$, and $A$ are matrices. Here, we will use the simplified version where $A$ is the tridiagonal stencil $[1,-2,1]$, i.e. it's the 2D discretization of the LaPlacian. The native code would be something along the lines of:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DiffEqBase.ODEProblem with uType Array{Float64,3} and tType Float64. In-place: true\n",
       "timespan: (0.0, 0.1)\n",
       "u0: [11.0168 11.0168 … 11.0168 11.0168; 11.0168 11.0168 … 11.0168 11.0168; … ; 11.0168 11.0168 … 11.0168 11.0168; 11.0168 11.0168 … 11.0168 11.0168]\n",
       "\n",
       "[12.1 12.1 … 12.1 12.1; 12.1 12.1 … 12.1 12.1; … ; 12.1 12.1 … 12.1 12.1; 12.1 12.1 … 12.1 12.1]"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Generate the constants\n",
    "p = (1.0,1.0,1.0,10.0,0.001,100.0) # a,α,ubar,β,D1,D2\n",
    "N = 100\n",
    "Ax = full(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))\n",
    "Ay = copy(Ax)\n",
    "Ax[2,1] = 2.0\n",
    "Ax[end-1,end] = 2.0\n",
    "Ay[1,2] = 2.0\n",
    "Ay[end,end-1] = 2.0                                \n",
    "\n",
    "function basic_version!(dr,r,p,t)\n",
    "  a,α,ubar,β,D1,D2 = p\n",
    "  u = r[:,:,1]\n",
    "  v = r[:,:,2]\n",
    "  Du = D1*(Ay*u + u*Ax)\n",
    "  Dv = D2*(Ay*v + v*Ax)\n",
    "  dr[:,:,1] = Du + a*u*u./v + ubar - α*u\n",
    "  dr[:,:,2] = Dv + a*u*u - β*v   \n",
    "end\n",
    "\n",
    "a,α,ubar,β,D1,D2 = p\n",
    "uss = (ubar+β)/α\n",
    "vss = (a/β)*uss^2\n",
    "r0 = zeros(100,100,2)\n",
    "r0[:,:,1] = uss+0.1*rand()\n",
    "r0[:,:,2] = vss\n",
    "\n",
    "prob = ODEProblem(basic_version!,r0,(0.0,0.1),p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this version we have encoded our initial condition to be a 3-dimensional array, with `u[:,:,1]` being the `A` part and `u[:,:,2]` being the `B` part. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  208.84 MiB\n",
       "  allocs estimate:  6753\n",
       "  --------------\n",
       "  minimum time:     17.634 s (3.46% GC)\n",
       "  median time:      17.634 s (3.46% GC)\n",
       "  mean time:        17.634 s (3.46% GC)\n",
       "  maximum time:     17.634 s (3.46% GC)\n",
       "  --------------\n",
       "  samples:          1\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark solve(prob,Tsit5())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While this version isn't very efficient,\n",
    "\n",
    "#### We recommend writing the \"high-level\" code first, and iteratively optimizing it!\n",
    "\n",
    "The first thing that we can do is get rid of the slicing allocations. The operation `r[:,:,1]` creates a temporary array instead of a \"view\", i.e. a pointer to the already existing memory. To make it a view, add `@view`. Note that we have to be careful with views because they point to the same memory, and thus changing a view changes the original values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = [0.220823, 0.296939, 0.979577, 0.850122]\n",
      "A = [0.220823, 2.0, 0.979577, 0.850122]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "4-element Array{Float64,1}:\n",
       " 0.220823\n",
       " 2.0     \n",
       " 0.979577\n",
       " 0.850122"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(4)\n",
    "@show A\n",
    "B = @view A[1:3]\n",
    "B[2] = 2\n",
    "@show A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that changing `B` changed `A`. This is something to be careful of, but at the same time we want to use this since we want to modify the output `dr`. Additionally, the last statement is a purely element-wise operation, and thus we can make use of broadcast fusion there. Let's rewrite `basic_version!` to ***avoid slicing allocations*** and to ***use broadcast fusion***:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  9.50 MiB\n",
       "  allocs estimate:  469\n",
       "  --------------\n",
       "  minimum time:     224.088 ms (0.00% GC)\n",
       "  median time:      498.966 ms (0.00% GC)\n",
       "  mean time:        532.433 ms (0.28% GC)\n",
       "  maximum time:     995.933 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          10\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function gm2!(dr,r,p,t)\n",
    "  a,α,ubar,β,D1,D2 = p\n",
    "  u = @view r[:,:,1]\n",
    "  v = @view r[:,:,2]\n",
    "  du = @view dr[:,:,1]\n",
    "  dv = @view dr[:,:,2]\n",
    "  Du = D1*(Ay*u + u*Ax)\n",
    "  Dv = D2*(Ay*v + v*Ax)\n",
    "  @. du = Du + a*u*u./v + ubar - α*u\n",
    "  @. dv = Dv + a*u*u - β*v   \n",
    "end\n",
    "prob = ODEProblem(gm2!,r0,(0.0,0.1),p)\n",
    "@benchmark solve(prob,Tsit5())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, most of the allocations are taking place in `Du = D1*(Ay*u + u*Ax)` since those operations are vectorized and not mutating. We should instead replace the matrix multiplications with `A_mul_B!`. When doing so, we will need to have cache variables to write into. This looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  4.01 MiB\n",
       "  allocs estimate:  397\n",
       "  --------------\n",
       "  minimum time:     101.338 ms (0.00% GC)\n",
       "  median time:      295.911 ms (0.00% GC)\n",
       "  mean time:        335.541 ms (1.94% GC)\n",
       "  maximum time:     898.436 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          15\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ayu = zeros(N,N)\n",
    "uAx = zeros(N,N)\n",
    "Du = zeros(N,N)\n",
    "Ayv = zeros(N,N)\n",
    "vAx = zeros(N,N)\n",
    "Dv = zeros(N,N)\n",
    "function gm3!(dr,r,p,t)\n",
    "  a,α,ubar,β,D1,D2 = p\n",
    "  u = @view r[:,:,1]\n",
    "  v = @view r[:,:,2]\n",
    "  du = @view dr[:,:,1]\n",
    "  dv = @view dr[:,:,2]\n",
    "  A_mul_B!(Ayu,Ay,u)\n",
    "  A_mul_B!(uAx,u,Ax)\n",
    "  A_mul_B!(Ayv,Ay,v)\n",
    "  A_mul_B!(vAx,v,Ax)\n",
    "  @. Du = D1*(Ayu + uAx)\n",
    "  @. Dv = D2*(Ayv + vAx)\n",
    "  @. du = Du + a*u*u./v + ubar - α*u\n",
    "  @. dv = Dv + a*u*u - β*v   \n",
    "end\n",
    "prob = ODEProblem(gm3!,r0,(0.0,0.1),p)\n",
    "@benchmark solve(prob,Tsit5())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But our temporary variables are global variables. We need to either declare the caches as `const` or localize them. We can localize them by adding them to the parameters, `p`. It's easier for the compiler to reason about local variables than global variables. ***Localizing variables helps to ensure type stability***."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  4.00 MiB\n",
       "  allocs estimate:  181\n",
       "  --------------\n",
       "  minimum time:     99.400 ms (0.00% GC)\n",
       "  median time:      282.105 ms (0.00% GC)\n",
       "  mean time:        357.870 ms (0.14% GC)\n",
       "  maximum time:     1.094 s (0.00% GC)\n",
       "  --------------\n",
       "  samples:          14\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = (1.0,1.0,1.0,10.0,0.001,100.0,Ayu,uAx,Du,Ayv,vAx,Dv) # a,α,ubar,β,D1,D2\n",
    "function gm4!(dr,r,p,t)\n",
    "  a,α,ubar,β,D1,D2,Ayu,uAx,Du,Ayv,vAx,Dv = p\n",
    "  u = @view r[:,:,1]\n",
    "  v = @view r[:,:,2]\n",
    "  du = @view dr[:,:,1]\n",
    "  dv = @view dr[:,:,2]\n",
    "  A_mul_B!(Ayu,Ay,u)\n",
    "  A_mul_B!(uAx,u,Ax)\n",
    "  A_mul_B!(Ayv,Ay,v)\n",
    "  A_mul_B!(vAx,v,Ax)\n",
    "  @. Du = D1*(Ayu + uAx)\n",
    "  @. Dv = D2*(Ayv + vAx)\n",
    "  @. du = Du + a*u*u./v + ubar - α*u\n",
    "  @. dv = Dv + a*u*u - β*v   \n",
    "end\n",
    "prob = ODEProblem(gm4!,r0,(0.0,0.1),p)\n",
    "@benchmark solve(prob,Tsit5())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could then use the BLAS `gemmv` to optimize the matrix multiplications some more, but instead let's devectorize the stencil. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  4.00 MiB\n",
       "  allocs estimate:  156\n",
       "  --------------\n",
       "  minimum time:     1.302 ms (0.00% GC)\n",
       "  median time:      1.605 ms (0.00% GC)\n",
       "  mean time:        2.088 ms (22.69% GC)\n",
       "  maximum time:     7.175 ms (62.19% GC)\n",
       "  --------------\n",
       "  samples:          2385\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = (1.0,1.0,1.0,10.0,0.001,100.0,N)\n",
    "function fast_gm!(du,u,p,t)\n",
    "  a,α,ubar,β,D1,D2,N = p\n",
    "\n",
    "  @inbounds for j in 2:N-1, i in 2:N-1\n",
    "    du[i,j,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +\n",
    "              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]\n",
    "  end\n",
    "\n",
    "  @inbounds for j in 2:N-1, i in 2:N-1\n",
    "    du[i,j,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +\n",
    "            a*u[i,j,1]^2 - β*u[i,j,2]\n",
    "  end\n",
    "\n",
    "  @inbounds for j in 2:N-1\n",
    "    i = 1\n",
    "    du[1,j,1] = D1*(2u[i+1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +\n",
    "            a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]\n",
    "  end\n",
    "  @inbounds for j in 2:N-1\n",
    "    i = 1\n",
    "    du[1,j,2] = D2*(2u[i+1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +\n",
    "            a*u[i,j,1]^2 - β*u[i,j,2]\n",
    "  end\n",
    "  @inbounds for j in 2:N-1\n",
    "    i = N\n",
    "    du[end,j,1] = D1*(2u[i-1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +\n",
    "           a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]\n",
    "  end\n",
    "  @inbounds for j in 2:N-1\n",
    "    i = N\n",
    "    du[end,j,2] = D2*(2u[i-1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +\n",
    "           a*u[i,j,1]^2 - β*u[i,j,2]\n",
    "  end\n",
    "\n",
    "  @inbounds for i in 2:N-1\n",
    "    j = 1\n",
    "    du[i,1,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +\n",
    "              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]\n",
    "  end\n",
    "  @inbounds for i in 2:N-1\n",
    "    j = 1\n",
    "    du[i,1,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +\n",
    "              a*u[i,j,1]^2 - β*u[i,j,2]\n",
    "  end\n",
    "  @inbounds for i in 2:N-1\n",
    "    j = N\n",
    "    du[i,end,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +\n",
    "             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]\n",
    "  end\n",
    "  @inbounds for i in 2:N-1\n",
    "    j = N\n",
    "    du[i,end,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +\n",
    "             a*u[i,j,1]^2 - β*u[i,j,2]\n",
    "  end\n",
    "\n",
    "  @inbounds begin\n",
    "    i = 1; j = 1\n",
    "    du[1,1,1] = D1*(2u[i+1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +\n",
    "              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]\n",
    "    du[1,1,2] = D2*(2u[i+1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +\n",
    "              a*u[i,j,1]^2 - β*u[i,j,2]\n",
    "\n",
    "    i = 1; j = N\n",
    "    du[1,N,1] = D1*(2u[i+1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +\n",
    "             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]\n",
    "    du[1,N,2] = D2*(2u[i+1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +\n",
    "             a*u[i,j,1]^2 - β*u[i,j,2]\n",
    "\n",
    "    i = N; j = 1\n",
    "    du[N,1,1] = D1*(2u[i-1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +\n",
    "             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]\n",
    "    du[N,1,2] = D2*(2u[i-1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +\n",
    "             a*u[i,j,1]^2 - β*u[i,j,2]\n",
    "\n",
    "    i = N; j = N\n",
    "    du[end,end,1] = D1*(2u[i-1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +\n",
    "             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]\n",
    "    du[end,end,2] = D2*(2u[i-1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +\n",
    "             a*u[i,j,1]^2 - β*u[i,j,2]\n",
    "   end\n",
    "end\n",
    "prob = ODEProblem(fast_gm!,r0,(0.0,0.1),p)\n",
    "@benchmark solve(prob,Tsit5())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, we can do other things like multithread the main loops, but these optimizations get the last 2x-3x out. The main optimizations which apply everywhere are the ones we just performed (though the last one only works if your matrix is a stencil. This is known as a matrix-free implementation of the PDE discretization).\n",
    "\n",
    "This gets us to about 20,000x faster than our original MATLAB/SciPy/R vectorized style code!\n",
    "\n",
    "The last thing to do is then ***optimize our algorithm choice***. We have been using `Tsit5()` as our test algorithm, but in reality this problem is a stiff PDE discretization and thus one recommendation is to use `CVODE_BDF()`. However, instead of using the default dense Jacobian, we should make use of the sparse Jacobian afforded by the problem. The Jacobian is the matrix $\\frac{df_i}{dr_j}$, where $r$ is read by the linear index (i.e. down columns). But since the $u$ variables depend on the $v$, the band size here is large, and thus this will not do well with a Banded Jacobian solver. Instead, we utilize sparse Jacobian algorithms. `CVODE_BDF` allows us to use a sparse Newton-Krylov solver by setting `linear_solver = :GMRES` (see [the solver documentation](http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#Sundials.jl-1), and thus we can solve this problem efficiently. Let's see how this scales as we increase the integration time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  33.31 MiB\n",
       "  allocs estimate:  588\n",
       "  --------------\n",
       "  minimum time:     25.352 ms (9.66% GC)\n",
       "  median time:      42.759 ms (44.63% GC)\n",
       "  mean time:        68.328 ms (62.93% GC)\n",
       "  maximum time:     181.947 ms (81.12% GC)\n",
       "  --------------\n",
       "  samples:          75\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prob = ODEProblem(fast_gm!,r0,(0.0,10.0),p)\n",
    "@benchmark solve(prob,Tsit5())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  11.48 MiB\n",
       "  allocs estimate:  1036\n",
       "  --------------\n",
       "  minimum time:     20.866 ms (0.00% GC)\n",
       "  median time:      24.317 ms (10.64% GC)\n",
       "  mean time:        24.123 ms (6.37% GC)\n",
       "  maximum time:     31.583 ms (17.76% GC)\n",
       "  --------------\n",
       "  samples:          208\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  282.49 MiB\n",
       "  allocs estimate:  4260\n",
       "  --------------\n",
       "  minimum time:     208.330 ms (0.00% GC)\n",
       "  median time:      450.052 ms (36.52% GC)\n",
       "  mean time:        475.101 ms (48.98% GC)\n",
       "  maximum time:     731.369 ms (70.71% GC)\n",
       "  --------------\n",
       "  samples:          12\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prob = ODEProblem(fast_gm!,r0,(0.0,100.0),p)\n",
    "@benchmark solve(prob,Tsit5())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  14.54 MiB\n",
       "  allocs estimate:  1282\n",
       "  --------------\n",
       "  minimum time:     26.585 ms (0.00% GC)\n",
       "  median time:      30.281 ms (8.68% GC)\n",
       "  mean time:        30.292 ms (6.16% GC)\n",
       "  maximum time:     37.065 ms (8.04% GC)\n",
       "  --------------\n",
       "  samples:          166\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What's happening is that, because the problem is stiff, the number of steps required by the explicit Runge-Kutta method grows rapidly, whereas `CVODE_BDF` is taking large steps. Additionally, the `GMRES` linear solver form is quite an efficient way to solve the implicit system in this case. This is problem-dependent, and in many cases using a Krylov method effectively requires a preconditioner, so you need to play around with testing other algorithms and linear solvers to find out what works best with your problem.\n",
    "\n",
    "## Conclusion\n",
    "\n",
    "Julia gives you the tools to optimize the solver \"all the way\", but you need to make use of it. The main thing to avoid is temporary allocations. For small systems, this is effectively done via static arrays. For large systems, this is done via in-place operations and cache arrays. Either way, the resulting solution can be immensely sped up over vectorized formulations by using these principles."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.2",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
